Bulletin of Mathematical Biology

Gliomas, a type of brain tumor, have become increasingly important in oncology as they are difficult to treat due to their location deep in the brain. While some research has been done on spatiotemporal prediction of future glioma growth--something that can aid in surgical resections of gliomas once a patient has been diagnosed--modeling efforts for pre-diagnostic gliomas remain limited. The lack of retrospective, pre-diagnostic data makes this a challenging task; yet, pre-symptomatic serum glucose levels in patients have been shown to have a relationship with the emergence of gliomas, motivating this area of research. In 2015, Sturrock et al. presented an ordinary differential equation model of pre-diagnostic glioma growth that describes glioma-glucose-immune interactions. This report reproduces the major findings of Sturrock et al., revising their model to incorporate more biological phenomena--namely vascularization and mutational burden--testing additional medically relevant patient scenarios, and providing an extended discussion on the implications of the model. In-silico simulations performed in this report provide further insight into models describing glioma-glucose-immune interactions, and how they can be expanded to incorporate physiologically relevant features. Future work is necessary to refine model parameters and validate predictions with the limited, albeit steadily growing, amount of longitudinal patient data.

Formulating tumor models that predict growth under therapy is vital for improving patient-specific treatment plans. In this context, we present our recent work on simulating non-small-scale cell lung cancer (NSCLC) in a simple, deterministic setting for two different patients receiving an immunotherapeutic treatment. At its core, our model consists of a Cahn-Hilliard-based phase-field model describing the evolution of proliferative and necrotic tumor cells. These are coupled to a simplified nutrient model that drives the growth of the proliferative cells and their decay into necrotic cells. The applied immunotherapy decreases the proliferative cell concentration. Here, we model the immunotherapeutic agent concentration in the entire lung over time by an ordinary differential equation (ODE). Finally, reaction terms provide a coupling between all these equations. By assuming spherical, symmetric tumor growth and constant nutrient inflow, we simplify this full 3D cancer simulation model to a reduced 1D model. We can then resort to patient data gathered from computed tomography (CT) scans over several years to calibrate our model. For the reduced 1D model, we show that our model can qualitatively describe observations during immunotherapy by fitting our model parameters to existing patient data. Our model covers cases in which the immunotherapy is successful and limits the tumor size, as well as cases predicting a sudden relapse, leading to exponential tumor growth. Finally, we move from the reduced model back to the full 3D cancer simulation in the lung tissue. Thereby, we show the predictive benefits a more detailed patient-specific simulation including spatial information could yield in the future. Author summaryLung cancer is one of the deadliest diseases, with low long-term survival rates. Its treatment is still very heuristic since patients respond to the same treatment plans differs significantly. Therefore, patient-specific models for predicting tumor growth and the treatment response are necessary for clinicians to make informed decisions about the patients therapy and avoid a trial and error based approach. We made a first small step in that direction by introducing a model for simulating cancer growth and its treatment inside a 3D lung geometry. In this model, we represented tumor cells by a volume fraction field that varies over space and time. We described their evolution by a system of partial differential equations, which include patient- and treatment-specific parameters capturing the different responses of patients to the therapies. Our simulation results were compared to clinical data and showed that we can quantitatively describe the tumors behavior with a suitable parameter set. This enabled us to change therapies in simulation runs and analyze how these changes could have impacted the patients health.

One of the most difficult challenges in cancer therapy is the emergence of drug resistance within tumors. Sometimes drug resistance can emerge as the result of mutations and Darwinian selection. However, recently another phenomenon has been discovered, in which tumor cells switch back and forth between drug-sensitive and pre-resistant states. Upon exposure to the drug, sensitive cells die off, and pre-resistant cells become locked in to a state of permanent drug resistance. In this paper, we explore the implications of this transient state switching for therapy scheduling. We propose a model to describe the phenomenon and estimate parameters from experimental melanoma data. We then compare the performance of continuous and alternating drug schedules, and use sensitivity analysis to explore how different conditions affect the efficacy of each schedule. We find that for our estimated parameters, a continuous therapy schedule is optimal. However we also find that an alternating schedule can be optimal for other, hypothetical parameter sets, depending on the difference in growth rate between pre-drug and post-drug cells, the delay between exposure to the drug and emergence of resistance, and the rate at which pre-resistant cells become resistant relative to the rate at which they switch back to the sensitive state.

Paclitaxel (taxol) is a commonly used chemotherapeutic that stabilizes microtubules to inhibit chromosome segregation during mitosis. Although it is used effectively to treat a wide variety of solid tumour types, it also impairs healthy cell proliferation leading to severe dose-limiting toxicities. Newer anti-mitotic drugs have been developed but these have so far failed to offer the same clinical benefit as paclitaxel, which begs the question of why this drug targets tumour cells so effectively? Here we develop a mathematical model of paclitaxel penetration and retention within 3D tumour environments following periodic drug-on/drug-off regimes typically used in the clinic. Our model suggests that during the drug-free periods, dense poorly-perfused tissue can retain paclitaxel for much longer than well-perfused tissue. This is due to paclitaxels ability to bind strongly to microtubules, which causes slower drug-release from densely packed tissue. Assuming that tumour cells are generally dense and less perfused than proliferative healthy tissues, this simple model suggests that tumour-selective killing could be achieved later in each chemotherapy cycle when the drug has otherwise cleared the healthy cell compartments. We use our model to optimize dosing regimens to allow paclitaxel to selectively kill tumour spheroids of different size, whilst sparing well-perfused healthy cells. Together, our model suggests paclitaxel could target key distinguishing features of many solid tumours: their size, 3D geometry and perfusion status. It is important to validate these predictions in cell models because, if correct, they could be harnessed to optimize paclitaxel use, to predict and enhance tumour responsiveness, and to develop newer drugs that are preferentially retained for longer within solid tumours.

Classical mathematical models of tumor growth have shaped our understanding of cancer and have broad practical implications for treatment scheduling and dosage. However, even the simplest textbook models have been barely validated in real world-data of human patients. In this study, we fitted a range of differential equation models to tumor volume measurements of patients undergoing chemotherapy or cancer immunotherapy for solid tumors. We used a large dataset of 1472 patients with three or more measurements per target lesion, of which 652 patients had six or more data points. We show that the early treatment response shows only moderate correlation with the final treatment response, demonstrating the need for nuanced models. We then perform a head-to-head comparison of six classical models which are widely used in the field: the Exponential, Logistic, Classic Bertalanffy, General Bertalanffy, Classic Gompertz and General Gompertz model. Several models provide a good fit to tumor volume measurements, with the Gompertz model providing the best balance between goodness of fit and number of parameters. Similarly, when fitting to early treatment data, the general Bertalanffy and Gompertz models yield the lowest mean absolute error to forecasted data, indicating that these models could potentially be effective at predicting treatment outcome. In summary, we provide a quantitative benchmark for classical textbook models and state-of-the art models of human tumor growth. We publicly release an anonymized version of our original data, providing the first benchmark set of human tumor growth data for evaluation of mathematical models. Author SummaryMathematical oncology uses quantitative models for prediction of tumor growth and treatment response. The theoretical foundation of mathematical oncology is provided by six classical mathematical models: the Exponential, Logistic, Classic Bertalanffy, General Bertalanffy, Classic Gompertz and General Gompertz model. These models have been introduced decades ago, have been used in thousands of scientific articles and are part of textbooks and curricula in mathematical oncology. However, these models have not been systematically tested in clinical data from actual patients. In this study, we have collected quantitative tumor volume measurements from thousands of patients in five large clinical trials of cancer immunotherapy. We use this dataset to systematically investigate how accurately mathematical models can describe tumor growth, showing that there are pronounced differences between models. In addition, we show that two of these models can predict tumor response to immunotherapy and chemotherapy at later time points when trained on early tumor growth dynamics. Thus, our article closes a conceptual gap in the literature and at the same time provides a simple tool to predict response to chemotherapy and immunotherapy on the level of individual patients.

Multiple Myeloma (MM) is a plasma cell cancer that occurs in the bone marrow. A leading treatment for MM is the monoclonal antibody Daratumumab, targeting the CD38 receptor, which is highly overexpressed in myeloma cells. In this work we model drug evasion via loss of CD38 expression, which is a proposed mechanism of resistance to Daratumumab treatment. We develop an ODE model that includes drug evasion via two mechanisms: a direct effect in which CD38 expression is lost without cell death in response to Daratumumab, and an indirect effect in which CD38 expression switches on and off in the cancer cells; myeloma cells that do not express CD38 have lower fitness but are shielded from the drug action. The model also incorporates competition with healthy cells, death of healthy cells due to off-target drug effects, and a Michaelis-Menten type immune response. Using optimal control theory, we study the effect of the drug evasion mechanisms and the off-target drug effect on the optimal treatment regime. We identify a general increase in treatment duration and costs, with varying patterns of response for the different controlling parameters. Several distinct optimal treatment regimes are identified within the parameter space. Author summaryIn this work we investigate a model of Multiple Myeloma, a cancer of the bone marrow, and its treatment with the drug Daratumumab. The model incorporates proposed mechanisms by which the cancer evades Daratumumab by reduced expression of the receptor CD38, which is the drug target and normally abundent in the cancer cells. The model includes an off-target effect, meaning that the drug treatment destroys some healthy cells alongside the targeted cancer cells. Both mechanisms can reasonably be expected to reduce the efficacy of the drug. We investigate the model using optimal control methods, which are used to find the drug dose over time which best balances the financial and health costs of treatment against cancer persistence, according to a specified cost function. We show that this drug resistence and off-target effect prolongs the optimal treatment and increase the burden of both the disease and drug. We analyse the distinct effects of the controlling parameters on each of these costs factors as well as the time course, and identify conditions under which extended treatment is required, with either intermittant treatment or a steady reduced dose. Extended treatment may be indefinite or for a fixed period.

Glucose is the principal metabolic fuel for the energy needs of most cell types. Upon infection, cytokines secreted by the immune system regulate redistribution of glucose to meet new metabolic needs associated with clearing the pathogen. We develop a mathematical model to describe the dynamics of such adaptation of metabolic pathways mediated by the immune response and its impact on the ability to clear pathogen and restore health. We find that cytokine-regulated redistribution of glucose resources in different tissues is critical for an effective immune response to pathogen as strictly clamping plasma glucose levels to homeostatic levels results in an ineffective immune response. By studying the effects of various parameters in our model, we describe how aberrant regulation of adaptation mechanisms affect outcomes of infection. Too high a glucose consumption rate by innate immune cells to mediate functions results in failure to clear pathogen. Pathogens with a very high replication rate can be controlled to low levels, but at a very high metabolic cost. Too low a pathogen replication rate allows the pathogen to hide from the immune system and rebound to high levels at later times. Finally, the strength of the innate immune response must be regulated to not be too high, not only to limit immunopathogenesis, but also for mediating an effective adaptive immune response.

Stochastic effects in cell growth and division drive variability in cellular volumes both at the single-cell level and at the level of growing cell populations. Here we consider a simple and tractable model in which cell volumes grow exponentially, cell division is symmetric, and its rate is volume-dependent. Consistently with previous observations, the model is shown to sustain oscillatory behaviour with alternating phases of slow and fast growth. Exact simulation algorithms and large-time asymptotics are developed and cross-validated for the single-cell and whole-population formulations of the model. The two formulations are shown to provide similar results during the phases of slow growth, but differ during the fast-growth phases. Specifically, the single-cell formulation systematically underestimates the proportion of small cells. More generally, our results suggest that measurable characteristics of cells may follow different distributions depending on whether a single-cell lineage or an entire population is considered.

Integrodifference equations (IDEs) are often used for discrete-time continuous-space models in mathematical biology. The model includes two stages: the reproduction stage, and the dispersal stage. The output of the model is the population density of a species for the next generation across the landscape, given the current population density. Most previous models for dispersal in a heterogeneous landscape approximate the landscape by a set of homogeneous patches, and allow for different demographic and dispersal rates within each patch. Some work has been done designing and analyzing models which also include a patch preference at the boundaries, which is commonly referred to as the degree of bias. Individuals dispersing across a patchy landscape can detect the changes in habitat at a neighborhood of a patch boundary, and as a result, they might change the direction of their movement if they are approaching a bad patch.\n\nIn our work, we derive a generalization of the classic Laplace kernel, which includes different dispersal rates in each patch as well as different degrees of bias at the patch boundaries. The simple Laplace kernel and the truncated Laplace kernel most often used in classical work appear as special cases of this general kernel. The form of this general kernel is the sum of two different terms: the classic truncated Laplace kernel within each patch, and a correction accounting for the bias at patch boundaries.

In this manuscript, we derive a closed form solution to the full Kermack and McKendrick integro-differential equations (Kermack and McKendrick 1927) which we call the KMES. We demonstrate the veracity of the KMES using independent data from the Covid 19 pandemic and derive many previously unknown and useful analytical expressions for characterizing and managing an epidemic. These include expressions for the viral load, the final size, the effective reproduction number, and the time to the peak in infections. The KMES can also be cast in the form of a step function response to the input of new infections; and that response is the time series of total infections. Since the publication of Kermack and McKendricks seminal paper (1927), thousands of authors have utilized the Susceptible, Infected, and Recovered (SIR) approximations; expressions putatively derived from the integro-differential equations to model epidemic dynamics. Implicit in the use of the SIR approximation are the beliefs that there is no closed form solution to the more complex integro-differential equations, that the approximation adequately reproduces the dynamics of the integro-differential equations, and that herd immunity always exists. However, the KMES demonstrates that the SIR models are not adequate representations of the integro-differential equations, and herd immunity is not guaranteed. We suggest that the KMES obsoletes the need for the SIR approximations; and provides a new level of understanding of epidemic dynamics.

Metabolism plays a complex role in the evolution of cancerous tumors, including inducing a multifaceted effect on the immune system to aid immune escape. Immune escape is, by definition, a collective phenomenon by requiring the presence of two cell types interacting in close proximity: tumor and immune. The microenvironmental context of these interactions is influenced by the dynamic process of blood vessel growth and remodelling, creating heterogeneous patches of well-vascularized tumor or acidic niches. We present a multiscale mathematical model that captures the phenotypic, vascular, microenvironmental, and spatial heterogeneity which shapes acid-mediated invasion and immune escape over a biologically-realistic time scale. We model immune escape mechanisms such as i) acid inactivation of immune cells, ii) competition for glucose, and iii) inhibitory immune checkpoint receptor expression (PD-L1) under anti-PD-L1 and sodium bicarbonate buffer therapies. To aid in understanding immune escape as a collective cellular phenomenon, we define immune escape in the context of six collective phenotypes (termed "meta-phenotypes"): Self-Acidify, Mooch Acid, PD-L1 Attack, Mooch PD-L1, Proliferate Fast, and Starve Glucose. Fomenting a stronger immune response leads to initial benefits but this advantage is offset by increased cell turnover that accelerates the emergence of aggressive phenotypes by inducing an evolutionary bottleneck. This model helps to untangle the key constraints on evolutionary costs and benefits of three key phenotypic axes on tumor invasion and treatment: acid-resistance, glycolysis, and PD-L1 expression. The benefits of concomitant anti-PD-L1 and buffer treatments is a promising treatment strategy to limit the adverse effects of immune escape. O_TEXTBOXSignificance statement In this work, we present a multi-scale mathematical model that captures the phenotypic, vascular, microenvironmental, and spatial heterogeneity which shapes acid-mediated invasion and immune escape over a biologically-realistic time scale. To aid in understanding immune escape as a collective cellular phenomenon, we introduce the concept of metaphenotypes: immune escape mechanisms that account for cellular phenotype and surrounding context. Metaphenotypes are defined by accounting for the following: cellular phenotype, microenvironmental factors, neighboring cell types, and immune-tumor interactions. These metaphenotypes provide insight into why targeting intratumoral pH with bicarbonate buffer is a synergistic combination treatment when paired with immune checkpoint blockade. C_TEXTBOX

Drug resistance in cancer is shaped not only by evolutionary processes but also by eco-evolutionary interactions between tumor subpopulations. These interactions can support the persistence of resistant cells even in the absence of treatment, undermining standard aggressive therapies and motivating drug holiday-based approaches that leverage ecological dynamics. A key challenge in implementing such strategies is efficiently identifying interaction between drug-sensitive and drug-resistant subpopulations. Evolutionary game theory provides a framework for characterizing these interactions. We investigate whether spatial patterns in single time-point images of cell populations can reveal the underlying game theoretic interactions between sensitive and resistant cells. To achieve this goal, we develop an agent-based model in which cell reproduction is governed by local game-theoretic interactions. We compute a suite of spatial statistics on single time-point images from the agent-based model under a range of games being played between cells. We quantify the informativeness of each spatial statistic and demonstrate that a simple machine learning model can classify the type of game being played. Our findings suggest that spatial structure contains sufficient information to infer ecological interactions. This work represents a step toward clinically viable tools for identifying cell-cell interactions in tumors, supporting the development of ecologically informed cancer therapies. Author summaryDrug resistance is a major challenge in cancer treatment, often leading to relapse despite initially successful therapy. While mutations are a key driver, ecological interactions between drug-sensitive and drug-resistant cells also play a critical role. These interactions are complex and dynamic, and few molecular biomarkers exist, making them difficult to study and account for in treatment planning. We use evolutionary game theory, a framework for quantifying interactions between cells, to investigate whether it is possible to infer these interactions using just a single time-point image of the cells. We develop an agent-based model where cells reproduce based on local interactions and quantify the resulting patterns in how cells are distributed across space using a suite of spatial statistics. We find that specific interaction types produce distinct spatial patterns that are evident in these metrics, and we train a simple machine learning model to classify the interaction type based on the metrics. Our results suggest that spatial data alone can offer valuable insights into tumor dynamics, potentially enabling more informed and adaptable cancer treatments based on eco-evolutionary principles.

Cytokine storm is a life-threatening inflammatory response that is characterized by hyperactivation of the immune system, and which can be caused by various therapies, auto-immune conditions, or pathogens, such as respiratory syndrome coronavirus 2 (SARS-CoV-2), which causes coronavirus disease COVID-19. While initial causes of cytokine storms can vary, late-stage clinical manifestations of cytokine storm converge and often overlap, and therefore a better understanding of how normal immune response turns pathological is warranted. Here we propose a theoretical framework, where cytokine storm phenomenology is captured using a conceptual mathematical model, where cytokines can both activate and regulate the immune system. We simulate normal immune response to infection, and through variation of system parameters identify conditions where, within the frameworks of this model, cytokine storm can arise. We demonstrate that cytokine storm is a transitional regime, and identify three main factors that must converge to result in storm-like dynamics, two of which represent individual-specific characteristics, thereby providing a possible explanation for why some people develop CRS, while others may not. We also discuss possible ecological insights into cytokine-immune interactions and provide mathematical analysis for the underlying regimes. We conclude with a discussion of how results of this analysis can be used in future research.

We assess the potential financial impact of future gene therapies by identifying the 109 late-stage gene therapy clinical trials currently underway, estimating the prevalence and incidence of their corresponding diseases, developing novel mathematical models of the increase in quality-adjusted life years for each approved gene therapy, and simulating the launch prices and the expected spending of these therapies over a 15-year time horizon. The results of our simulation suggest that an expected total of 1.09 million patients will be treated by gene therapy from January 2020 to December 2034. The expected peak annual spending on these therapies is $25.3 billion, and the total spending from January 2020 to December 2034 is $306 billion. We decompose their annual estimated spending by treated age group as a proxy for U.S. insurance type, and consider the tradeoffs of various methods of payment for these therapies to ensure patient access to their expected benefits.

We extend Bayesian Logistic Regression to model the dose-toxicity relationship in the setting of phase I dose-escalation/ dose-finding trials for cancer immunotherapies. Immunotherapy drugs are associated with Cytokine Release Syn-drome, a systemic immune system reaction that can be mitigated when initial lower doses of the drug are administered to generate immune tolerance. This changes the classic dose-finding problem of determining an optimal safe dose, to a more complex problem where the search is for both the optimal safe dose and the dose regimen that allows patients to quickly and safely reach that dose without CRS. As part of solving this methodological challenge, we show how to jointly model CRS and non-CRS toxicities, which have distinct mechanisms, while controlling for the overall toxicity rate to make dose-escalation decisions.

The Fisher-KPP model supports travelling wave solutions that are successfully used to model numerous invasive phenomena with applications in biology, ecology, and combustion theory. However, there are certain phenomena that the Fisher-KPP model cannot replicate, such as the extinction of invasive populations. The Fisher-Stefan model is an adaptation of the Fisher-KPP model to include a moving boundary whose evolution is governed by a Stefan condition. The Fisher-Stefan model also supports travelling wave solutions; however, a key additional feature of the Fisher-Stefan model is that it is able to simulate population extinction, giving rise to a spreading-extinction dichotomy. In this work, we revisit travelling wave solutions of the Fisher-KPP model and show that these results provide new insight into travelling wave solutions of the Fisher-Stefan model and the spreading-extinction dichotomy. Using a combination of phase plane analysis, perturbation analysis and linearisation, we establish a concrete relationship between travelling wave solutions of the Fisher-Stefan model and often-neglected travelling wave solutions of the Fisher-KPP model. Furthermore, we give closed-form approximate expressions for the shape of the travelling wave solutions of the Fisher-Stefan model in the limit of slow travelling wave speeds, c << 1.

Understanding and predicting the eco-evolutionary dynamics of cancer requires identifying mathematical models that best capture tumor growth and treatment response. In this study, we fit a family of two-population models to in-vitro data from non-small cell lung cancer (NSCLC), tracking drug-sensitive and drug-resistant cells under varying environmental conditions. The dataset, originally presented by Kaznatcheev et al., includes conditions with and without the drug Alectinib and cancer-associated fibroblasts (CAFs). We compare combinations of growth models (logistic, Gompertz, von Bertalanffy) and drug efficacy terms (Norton-Simon, linear, ratio-dependent) to identify which best explains the observed dynamics. Our models incorporate density dependence, frequency-dependent competition, and drug response, enabling mechanistic interpretation of tumor cell interactions. The logistic model with ratio-dependent drug efficacy best fits monoculture data. Using growth parameters from monocultures, we infer inter-type competition coefficients in co-cultures. We find that growth rate and carrying capacity are stable across CAF conditions, while competition and drug efficacy parameters shift, altering interaction dynamics. Notably, CAFs promote coexistence between resistant and sensitive cells, whereas Alectinib induces competitive exclusion. Our results underscore the need to evaluate both model fit and biological plausibility to guide therapeutic modeling of cancer. Author summaryHow cancer cells grow, compete, and respond to treatment depends not only on the drug, but also on their ecological context, including interactions with other cells and components of the tumor microenvironment. In this study, we explore how different mathematical models capture the behavior of non-small cell lung cancer (NSCLC) cells under various conditions. We focus on two cell populations: one sensitive to treatment, and one resistant. Using in-vitro data, we compare growth models and drug response types to identify which model best explains the observed population dynamics. We also investigate how different factors, such as the drug Alectinib and cancer-associated fibroblasts (CAFs), change the way cancer cells interact. Our game-theoretic approach allows us to quantify how these external conditions affect competition between cell types, revealing when resistant and sensitive cells can coexist. These findings contribute to a deeper understanding of tumor ecology and may support the development of adaptive cancer therapies that anticipate evolutionary responses to treatment.

In pediatric acute lymphoblastic leukemia (ALL), central nervous system (CNS) involvement is staged using white blood cell (WBC) counts from lumbar punctures (LPs), microscopy of LP-derived cells, and CNS imaging. CNS staging informs the need for additional intrathecal chemotherapy, which can result in side-effects including significant neurotoxicity. However, nearly 20% of LPs are traumatic, or contaminated by peripheral blood. In these cases, the Steinherz-Bleyer (S-B) algorithm is used instead of WBC counts to determine CNS involvement. The intuition for this algorithm has not been presented. In this work, we conceptualize LP lab values as a mixture of blood and cerebrospinal fluid (CSF) and present this mixture as a convex combination problem. Then, we demonstrate that the S-B algorithm asks whether the CSF to blood WBC ratio is at or above that of the mixing ratio (i.e., contaminated to uncontaminated ratio). Additionally, we derive an expression for estimating the true but unobserved CSF WBC in traumatic LP cases such that the existing atraumtic CSF WBC guidelines may be used. Finally, we present a Bayesian approach to incorporate non-zero CSF red blood cell (RBC) counts and suggest that this biologically-motivated assumption underlying S-B is likely not clinically relevant for the majority of patients.

Glioblastoma is a highly aggressive and difficult-to-treat brain cancer that resists conventional therapies. Recent advances in chimeric antigen receptor (CAR) T-cell therapy have shown promising potential for treating glioblastoma; however, achieving optimal efficacy remains challenging due to tumor antigen heterogeneity, the tumor microenvironment, and T-cell exhaustion. In this study, we developed a mathematical model of CAR T-cell therapy for glioblastoma to explore combinations of CAR T-cell treatments that take into account the spatial heterogeneity of antigen expression. Our hybrid model, created using the multicellular modeling platform PhysiCell, couples partial differential equations that describe the tumor microenvironment with agent-based models for glioblastoma and CAR T-cells. The model captures cell-to-cell interactions between the glioblastoma cells and CAR T-cells throughout treatment, focusing on three target antigens: IL-13R2, HER2, and EGFR. We analyze tumor antigen expression heterogeneity informed by expression patterns identified from human tissues and investigate patient-specific combination CAR T-cell treatment strategies. Our model demonstrates that an early intervention is the most effective approach, especially in glioblastoma tumors characterized by mixed antigen expression. However, in tissues with clustered antigen patterns, we find that sequential administration with specific CAR T-cell types can achieve efficacy comparable to simultaneous administration. In addition, spatially targeted delivery of CAR T-cells to specific tumor regions with matching antigen is an effective strategy as well. Our model provides a valuable platform for developing patient-specific CAR T-cell treatment plans with the potential to optimize scheduling and locations of CAR T-cell injections based on individual antigen expression profiles.

Drug resistance is a major challenge for cancer therapy. While resistance mutations are often the focus of investigation, non-genetic resistance mechanisms are also important. One such mechanism is the presence of relatively high fractions of cancer stem cells (CSCs), which have reduced susceptibility to chemotherapy, radiation, and targeted treatments compared to more differentiated cells. The reasons for high CSC fractions (CSC enrichment) are not well understood. Previous experimental and mathematical modeling work identified a particular feedback loop in tumors that can promote CSC enrichment. Here, we use mathematical models of hierarchically structured cell populations to build on this work and to provide a comprehensive analysis of how different feedback regulatory processes that might partially operate in tumors can influence the stem cell fractions during somatic evolution of healthy tissue or during tumor growth. We find that depending on the particular feedback loops that are present, CSC fractions can increase or decrease. We define characteristics of the feedback mechanisms that are required for CSC enrichment to occur, and show how the magnitude of enrichment is determined by parameters. In particular, enrichment requires a reduction in division rates or an increase in death rates with higher population sizes, and the feedback mediators that achieve this can be secreted by either CSCs or by more differentiated cells. The extent of enrichment is determined by the death rate of CSCs, the probability of CSC self-renewal, and by the strength of feedback on cell divisions. Defining these characteristics can guide experimental approaches that aim to screen for and identify feedback mediators that can promote CSC enrichment in specific cancers, which in turn can help understand and overcome the phenomenon of CSC-based therapy resistance.